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We present a simple analytical approach to study anharmonic effects in single layer, bilayer, and 
multilayer graphene. The coupling between in plane and out of plane modes leads to negative 
Griineisen coefficients and negative thermal expansion. The value of the thermal expansion coeffi- 
cient depends on the coupling to the substrate. The bending rigidity in bilayer graphene shows a 
crossover between a long wavelength regime where its value is determined by the in plane elastic 
properties and a short wavelength regime where its value approaches twice that of a single layer. 



Introduction. 

To realize the full potential of graphene layers in 
promising applications, like the design of fast electronic 
devices or sensitive and accurate molecular detectors, it 
is important to reach a thorough understanding of its 
properties down to the atomic level. [HE] At T= K, in 
the absence of defects, the carbon bond on the graphene 
layer is well understood in terms of the formation of three 
in-plane localized strong sp 2 bonds, and a fourth delocal- 
ized, out-of-plane, 7r-like bond. [<2 13 The optimum geo- 
metrical configuration is achieved by a honeycomb lattice 
formed by two equivalent sublattices displaying P6/mmm 
symmetry. The corresponding electronic structure shows 
bands dispersing linearly around the Fermi energy that 
are responsible for the fast and efficient transport of car- 
riers. Both experimentally and theoretically, [SHZj it is 
shown that this kind of arrangement results in a mate- 
rial with the largest in-plane elastic constants known yet. 
Therefore the 2D perfect flat layer makes the most sta- 
ble configuration since deviations from a common plane 
requires a significant amount of energy. Any departure 
from such a scenario affects greatly the atomic scale prop- 
erties of the layer and must be understood in order to 
efficiently exploit graphene's properties. Different rea- 
sons, however, might be invoked for a two-dimensional 
graphene layer to adopt a certain corrugation at differ- 
ent scales. First, at a non-zero temperature a thermo- 
dynamic argument implies the impossibility for a perfect 
2D layer to exist in 3D.[H I5HT2"] Second, defects like ad- 
sorbed impurities, vacancies, etc, create local corruga- 
tions at the atomic scale[13 that propagate via the elas- 
tic properties of the lattice originating long-range corre- 
lations. Finally, external applied stresses related to con- 
ditions on the boundary make graphene to bend and to 
corrugate; an interesting point to study since the growth 
of graphene layers on different supporting substrates im- 
plies mismatches that introduce all kind of stresses that 
have been observed to originate a highly complex and 
corrugated landscape. [TH [15] In this work, we analyze 
a simple model based on the theory of elasticity to ob- 



tain physical insight on the Griineisen coefficients and 
the thermal expansion coefficient of graphene. [TBI which 
can be compared to atomistic models based on ab initio 
Density Functional Theory that yields a realistic quan- 
titative description of bending modes and corrugations 
appearing at the atomic scale |17j. 
Single layer graphene. 

We study anharmonic effects using the continuum the- 
ory of elasticity. We extend previous analyses [HI [19], us- 
ing the standard theory of free standing membranes [U03- 
ITT] . The Hamiltonian is [2D] 
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where p is the mass density, u is the two dimensional 
displacement vector, h is the displacement in the out of 
plane direction, k is the bending rigidity, and A and p 
are elastic Lame coefficients. For graphene, we have [7] 
k « 1 cV, A = 2 eV A" 2 and p = 10 eV A" 2 . 

We study the modes associated to the out of plane 
displacements. If we assume that there are no in plane 
tensions, diUj = 0, and we neglect the quartic terms in 



h, we obtain = y n |q| /p. This is the well known 
dispersion relation for out of plane flexural modes. We 
now analyze how these frequencies are modified when the 
in plane lattice constant is modified. An isotropic change 
of the lattice constant by a factor u can be included in 
the Hamiltonian, eq. [TJ by assuming that d x u x = d y u y = 



2 



u. The effective Hamiltonian for h, expanded to second 
order, becomes 
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The new frequencies of the ficxural phonons are 
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The derivative of the phonon frequency with respect to 
a change in the area of the unit cell A is 



7q 



A duj^ 
Wq- dA 



1 duj^ 



2wej du 



a=0 



x + n 

'2«|q| 2 



(4) 



where 73 is the Griineisen parameter. We obtain a neg- 
ative Griineisen parameter for all low frequency flexural 
modes, which diverge for |q| -> as |q| . This expres- 
sion is valid for momenta much smaller than the inverse of 
the interatomic spacing a, |q| <C a~ 1 . This result is con- 
sistent with a number of numerical calculations, which 
show negative Griineisen coefficients for flexural modes, 
which tend to diverge at low momenta[6l 171 |2T| 122] . 

Within the harmonic approximation, the estimate of 
the Griineisen parameters in eq. [4] allows us to obtain 
the thermal expansion coefficient [I] 
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where A is the area of the unit cell and we take into ac- 
count that the two-dimensional bulk modulus B = A+jU. 
The sum (in thermodynamic limit is replaced by an inte- 
gral) in the right-hand side of eq. ^ is divergent at small 
q which is the consequence of inapplicability of the har- 
monic approximation at small q where renormalization 
of effective bending rigidity and elastic modulii become 
relevant. The crossover wave vector is [H[5] 
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where Y = 4y«(A + /Lt)/(A + 2fi) is the two-dimensional 
Young modulus. Note that the corresponding phonon 
frequency lies deeply in the classical region: 
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where m and M are electron mass and mass of carbon 
atom, respectively. With the logarithmic accuracy, 
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where qr is the thermal wave vector satisfying the con- 
dition hu}(qr) = k B T. From the estimation in eq. pb, 
we obtain a w — K _1 , a quite good estimation for 
so oversimplified model (cf. Refs. [7|). Here we as- 
sume that the temperature is smaller than the maximal 
energy of the flexural phonon, T m « 15 THz w 700 K [B], 
otherwise one needs to add the factor T m /T under the 
argument of logarithm in eq. Q. 

Due to eq. ([7]) phonons relevant for the thermal ex- 
pansion coefficient can be considered as classical at any 
temperatures. This allows us to repeat the calculation 
of a taking into account anharmonic effects. Due to eq. 
([2]) and Hellmann-Feynman theorem the derivative of the 
free energy J- with respect to the deformation at u = 
can be rigorously expressed via the correlation function 
of out-of-plane displacements: 
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and via the anharmonic self energy S (~q*): 
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The latter can be estimated from the condition that at 
q = q* both terms in the denominator in eq. ( 10 1 are of 
the same order of magnitude [2"5] : 
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where 77 rs 0.85 is the exponent of renormalization of the 
bending rigidity; the numerical factor A was calculated 
within the self-consistent screening approximation [24] ; 
it was also shown that this approximation agrees quite 
well with the atomistic Monte Carlo simulations. Substi- 
tuting eq. (Ill into eq. (10 1 and further into eq. Q one 



can calculate the thermal expansion coefficient 
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with anharmonic effects taken into account. With the 
logarithmic accuracy, the result coincides with eq. {[sj) . 
Thus, the contribution of flexural mode to the thermal 
expansion coefficient is always negative and tempera- 
ture independent up to T ~ T m 700 K; at higher 
temperatures it depends on the temperature logarithmi- 
cally. This means that the inversion of sign of the ther- 
mal expansion coefficient at high temperature found in 
atomistic simulations [7] is due to contributions of other 
phonon modes. 

This justifies the use of quasiharmonic approximation 
to estimate the contribution of flexural phonons to the 
thermal expansion. Further we will consider only this 
approximation. 

Finally, from eq. [3] for the phonon frequencies we 
can estimate the momentum q c for which the value of 
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uri becomes negative for negative u. We obtain q c — 
y/[(X + h)\u\]/k. For u = -0.04 we find q c 0.6 A" 1 . 
Graphene on a substrate. 

The flcxural modes of graphene on a substrate are 
modified by the coupling to the substrate. The leading 
effect at long wavelengths can be analyzed by consid- 
ering the interaction energy per unit area between the 
graphene layer and the substrate, V su b s (h su b s ), where 
h su b s is the distance to the substrate. The disper- 
sion relation for the flexural modes becomes u\m ~ 

^/{^\ i + V"{h eq )]/p = ^{cj^+lj 2 , where h eq is the 

equilibrium distance. We can get an estimate of V"(h eq ) 

V"{h eq ) = (13) 

where V(h eq ) is the binding energy per unit area of 
graphene to the substrate, and do is a length scale such 
that do < h eq . The binding energy between graphene and 
the substrate depends on the precise attraction mech- 
anism |25| between the two materials, and it is likely 



bound by the van der Waals interactions. A reasonable 
range of values is 5 — 50 meV A~ 2 . For do ~ 2 A, we find 
hu>o ~ 1 — 4 meV. The value of ojq provides a cutoff in 
the expression for the thermal expansion, eq. [5] Hence, 
the negative contribution of the flcxural modes is reduced 
at temperatures such that T ps (hu )/kB ~ 10 — 40 K. 
For ksT Hujo the thermal expansion of graphene on 
a substrate should be similar to that of free standing 
graphene. At room temperature and V(h eq ) < 50 meV 
A~ 2 , the anharmonic momentum cutoff, q* (see eq. [fj]) is 
such that q* <C \V{h eq )/ Kd^] 1 ^, and the thermal expan- 
sion of graphene should be independent of the substrate. 

Bilayer graphene. 

In a discrete stack made of weakly coupled slabs we 
can expand the interlayer coupling assuming that the dis- 
placements vary slowly as a function of the two dimen- 
sional coordinate r, u 2 z — > (u nz — u n j r \ z /d) 2 ,u 2 , z +u 2 z — ► 

|(r n+ i — r„_i)/2 + V||U„ Z | , where d is the distance be- 
tween the layers. 

For two layers, an approximate expression is 
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For an infinite three dimensional stack, the parameters 
<7i, g 2 and g% define a continuum model like the one in 
eq. ([20) with c 33 = gi/d, c i3 = g 3 /d and c 44 = gild. 

If we assume that g 2 = 0, the in plane and out of plane 
modes are decoupled. The equations of motion for the 
out of plane modes are 

p 2D d 2 t u lz = -k (d 2 xx + d 2 yy ) 2 u lz - (u lz - u 2z ) 
P2Dd 2 t u 2z = k (d 2 xx + d 2 yy f u 2z - -± {u 2z - u 2z ) (15) 



where p is the two dimensional mass density. In mo- 
mentum space, we obtain two flexural modes, w+(k) = 
y^/pfc 2 , w_(k) = ^J{nk i + 2 gi /d 2 )/p, where k = |k|. 



For the phonon frequencies are obtained from 

the diagonalization of the 6x6 matrix. It can be split 
into two 3x3 matrices by using the combinations Yi — 
ir*2,Mi z = ^u 2z . The low energy modes are given by 
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The out of plane displacement couples to the longi- tudinal acoustical phonons. At low momenta we have 
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The quartic term in this expression is consistent with 
the continuum analysis described below. The flexural 
modes acquire a contribution which is independent of 
the parameter 172, and which scales with the three di- 
mensional bulk modulus and with d 3 , as the relation 
between two and three dimensional Lame coefficients is 
A 2 d,M2D oc Ad, [id. 

For graphene, k <C (A 2 d + 2/i 2 D)d 2 , so that the second 
term dominates in eq. (17 1. The bending rigidity of a 



bilayer should be significantly larger than that of a single 
layer, provided that the interlayer shear rigidity g 2 =/= 0. 
Using again A = 2 eV A" 2 and fi = 10 eV A" 2 , g 2 = 0.03 
eV and d = 3.3 A, we find a crossover from a high to a 
low value of the flexural rigidity at a length t = fc _1 rj 
55 A. Note that the atomistic simulations for finite-size 
crystallites in Ref. |26]deal with a larger k region giving 
approximately the same values for the bending rigidity 
(per layer) for single-layer and bilayer graphene. The 
dispersion of the flexural phonons is shown in Fig. [T] 

Griineisen coefficients in a bilayer. 

By applying an in plane strain, u, the frequencies in 



eq. (17 1 are reduced by 
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This expression gives a Griineisen parameter 

A 2 £> + ^2D A 2 d + M2-D 



Ik 



2 [k + (A 2 d + 2p 2D )d?] fc 2 2(X 2D + 2p 2D )d?k 2 

(19) 



This value is lower than the corresponding expression 
for single layer graphene, so that the negative expansion 
coefficient is reduced in a graphene bilayer. The analy- 
sis probably can be extended to graphite, although the 
dispersion of the out of plane modes will no longer be 
quadratic. 

Multilayered graphene. Continuum model. 
The elastic energy of a slab which is isotropic in the 
x — y plane can be written as: 




FIG. 1: Log-log plot of the dispersion of the flexural modes 
in a graphene bilayer using the parameters described in the 
text. The two straight lines correspond to the long and short 
wavelength limits discussed in the text. 



where we use te notation cy for the elastic constants 
instead of Lame coefficients. 

We assume that the slab is sufficiently narrow so that 
the stresses at the top and bottom surface do not differ 
much. The boundary conditions are[20j 
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From these equations, we obtain 
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Finally, the frequencies of the flexural modes are given 



by 
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This expression does not depend on the value of C44, but 
the value of this parameter must be different from zero, 



in order for eq. (21) to be valid, in agreement with the 



analysis carried out earlier for the bilayer. 
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